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ABSTRACT 

If the two planets in the HAT-F-IS system are coplanar, the orbital states provide a probe of the 
internal planetary structure. Frevious analyses of radial velocity and transit timing data of the system 
suggested that the observational constraints on the orbital states were rather small. We reanalyze the 
available data, treating the jitter as an unknown MCMC parameter, and find that a wide range of jitter 
values are plausible, hence the system parameters are less well constrained than previously suggested. 
For slightly increased levels of jitter (~ 4.5 ms^^) the eccentricity of the inner planet can be in the 
range < tinner < 0.07, the period and eccentricity of the outer planet can be 440 < Pouter < 470 
days and 0.55 < Couter < 0.85 respectively, while the relative pericenter alignment, ry, of the planets 
can take essentially any value —180° < rj < +180°. It is therefore difficult to determine whether 
Sinner and T] have evolved to a fixed-point state or a limit cycle, or to use Cmner to probe the internal 
planetary structure. We perform various transit timing variation (TTV) analyses, demonstrating 
that current constraints merely restrict Couter < 0.85, and rule out relative planetary inclinations 
within ~ 2° of Vei — 90°, but that future observations could significantly tighten the restriction on 
both these parameters. We demonstrate that TTV profiles can readily distinguish the theoretically 
favored inclinations of irei = 0°&45°, provided that sufficiently precise and frequent transit timing 
observations of HAT-F-13b can be made close to the pericenter passage of HAT-F-lSc. We note the 
relatively high probability that HAT-F-13c transits and suggest observational dates and strategies. 
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1. INTRODUCTION 

The HAT-P-13 system ()Bakos et al.|[20Q9l ) was the first extrasolar planetary system to be discovered in which both 
a transiting planet and an additional confirmed companion were known to coexist. Since this initial discovery, further 
transit-plus-companion systems have been discovered: Co RoT-7 (Queloz et al. 2 0C)9[) and HAT-P -7 (Narita et al. 201(J) 
as well as the recent multi-transit systems from Kepler ([Steffen et al.ll2010l : IHolman et al.ll2010 l ). while a number of 
other transiting systems display RV trends symptomatic of outer companions (E.g. HAT-P-11. iBakos et al.ll2010D . 

Such systems are of interest for a number of reasons. The first relates to the observable effects which arise from 
interactions between the two planets. The gravitational interaction between multiple planets causes the planetary 
orbits to be perturbed away from Keplerian ellipses. When one of the planets is transiting, these perturbations mean 
that the duration of, and the period between, successive transits will not be strictly constant. It has been calculated 
that observations of such transit timing variations (TTVs) and transit duration variations (TPy s ) would allow th e 
(inferred) dete ction of terrestrial-mass planets in Hot- Jupiter systems (IHolman and Murravl 120051 lAgol et al.l [20051) . 
trojan planets ()Ford and Holma"nll2007D and exoplanet moons IKipDind (|2009dia ). 

In addition, accur ate measure ments of transiting systems can allow us to observationally determine a huge range of 
system parameters (jWin n"2009') . one such parameter being the tidal Love number of the planet (if the system archi- 
tec ture is conveii i ent - | Wu and Goldreich (2002); Mardhng (2007); Ragozzine and Woh (2009)). In an investigation 
bv lBatvgin et al.l (12009) , it is demonstrated for the coplanar case that the HAT-P-13 system is indeed such a system, 
and a relationship is found between the Love number and the eccentricity of the inner planet. 

However, the determination of the Love number depends on the assumption that the inner planet tends towards a 
quasi-fixed point in (cinnenV) space, where emner is the eccentricity of the inner planet and r] = w outer — v^i nne.r is the 
differe nce in the alignment of the longitude of pericenters of the outer and inner planets. The recent work bv lMardlingl 
(|2010( ) looked at the evolution of a general non-coplanar, two-planet system, in which the angular momentum of the 
outer planet dominates the angular momentum budget of the system, and revealed that in such a system the inner 
planet does not tend to a fixed point, but instead tends to a limit cycle, with Cmner and t] constantly sampling along 
a closed trajectory. 

The approach to the limit cycle is found to be strongly dependent on the relative inclination between the two planets, 
ireV- The average value of Cinner around the cycle decreases and the limit cycle amplitude increases with increasing 
irei\ Limit cycle behaviour only exists for 0° < Vei < 33°, 46° < Vei < 54°, 126° < Vei < 134° and 147° < Vei < 180° 
For the regions 33° < irei < 46° and 134° < irei < 147°, 77 circulates and no limit cycle exists. For the region 
54° < irei < 126°, the effects of Kozai oscillations combined with tidal dissipation act to move the relative inclinations 
back towards irei — 54 or irei = 126 for prograde or retrograde orbits respectively, meaning that the system simply 
cannot exist with 54° < irei < 126° for a tidal dissipati on factor, Q < 10^. For the range of possible emner and 
Rinner (the radius of the inner planet) reported by Batvgi n et al.l (|2009t l. certain relative inclinations can be ruled out, 
suggesting that (for the prograde case), the system either has i < 10°, or i 45°. 

It is apparent from the analyses in both iMardlin^ ()2010[ ) and iBatvgin et al.l ()2009D that the conclusions one can 
derive regarding the possible state of the system rather sensitively depend on (i) the measured eccentricity of the inner 
planet, (ii) the relative pericenter of the two planets and (iii) the essentially unknown relative inclination between the 
two planets. It is our aim in this paper to try and understand in more detail what the current observational constraints 
on these quantities are. 

We start by re-evaluating the published system parameters for HA T-P-13 from IBakos et al.l ()2009t ) (henceforth B09) 
as well as those from the expanded analysis bv I Winn et al.l (|2010[ ) (henceforth WIO), performing a Markov chain 
Monte-Carlo (MCMC) investigation that concentrates on the effect of assumptions about jitter and how to include 
it within a MCMC analysis. We use jitter as a mathematical parameter to quantify the magnitude of unmodeled 
variations in the radial velocity observations. Jitter can be due to undetected planets, stellar activity or unrecognized 
noise in the instrument or data analysis pipeline. In the case of a complete model for the planetary system and ideal 
measurements, the jitter would reduce to the "stellar jitter" . Stellar jitter is the noise introduced into radial velocity 
(RV) measurements by unknown changes in the surface of the observed star, driven by sun spots, bulk flows and other 
inhomogeneities on the stellar surface. Several investigation s have examined this phenomenon, trying t o correlate th e 
magnitude of the jitter with observable stellar parameters (jSaar and Donahuill997l: ISaaFeTall Il998l: TWrightllMl . 
They find that stars of a given type can have a wi de range of jitter values, with stars similar to HAT-P-13 having 
jitters in the range 2 — 15 ms^^ fsee iWrighO ()2005[ ) and section IXD for further details). 

In the papers of B09 and WIO, subsequent to the MCMC analysis of the fitted planetary parameters, jitter levels of 
cTj — 3.0 m and aj = 3.4 m respectively were required in order to give reduced values of 1 (see ^for further 
discussion). We wish to understand the effects of changing this methodology, and in particular to find out to what 
degree including the jitter as an MCMC model parameter loosens the constraints on quantities such as the eccentricity 
of the inner planet and the pericenter alignment of the two planets. 

In a rigorous Bayesian analysis, the jitter should be treated as a model parameter simultaneously with the mass 
and orbital parameters. In this paper we show that correlations between the jitter and orbital parameters can lead 
to erroneous inferences if the jitter is held fixed during the modeling process. Thus, one should include the jitter in 
MCMC analyses alongside all of the other parameters that one is attempting to model. In this manner, one can arrive 
at a consistent statistical interpretation of ones knowledge of the observables in the system. In a Bayesian analysis, 
one must explicitly state assumptions for the prior distributions. While un- modeled observations (e.g., stellar color 
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or temperature) can influence the choice of prior (e.g., F stars are more likely to have a jitter exceeding 10 m/s), the 
choice of prior for the jitter must not be influenced by the radial velocities themselves. E.g., choosing a prior for the 
jitter based on the residuals to a fit results in double-counting the radial velocity data and can result in misleadingly 
small uncertainties. 

Further to this MCMC investigation, we go on to add an analysis of the TTVs for HAT-P-13. We do this to 
try and ascertain whether (i) combining the MCMC analysis with the TTV analysis can further restrict the range 
of parameter space available to the observed system quantities (eccentricities, alignments, etc), and (ii) we wish to 
understand whether the TTVs can provide some insight into the relative inclination of the planets (Relative planetary 
inclinations have previo usly been shown to be important in determining the expected TTVs in some svstems. iNesvorir^ 
l2009t IPavneind1[20Tnl ). 

We proceed in this paper as follows: In section [2] we outline the numerical methods we use to conduct our MCMC 
and TTV analyses; In section [3] we look at the effects of jitter in an MCMC analysis of the orbital elements of the 
planets in the HAT-P-13 system, and combine this with TTV constraints to understand whether this allows a more 
nuanced determination of the system parameters; In section |4] we look more generally at the TTVs in the HAT-P-13 
system, focusing in detail on the effects that relative planetary inclinations may have on the expected TTVs; After 
this, we move on in section[5]to consider the potential for future observations (both transit and RV) to better constrain 
the planetary orbits; and finally, in section [6] we present a summary and discussion of our conclusions. 

2. METHODOLOGY 
2.1. Radial Velocity & Transit Observation Analysis 

We analyze the radial velocity and transit observations using a Bayesian framework following iFordI ()2005l I2006D . 
We assume priors that are uniform in log of orbital period, eccentricity, argument of pericenter, mean anomaly at 
epoch, and the velocity zero-point. For the velocity amplitude (K) and jitter {(Jj), we adopt a prior of the form 
p{x) = (x + Xn)~^\loq(l + x/x„)]~ \ with Ko = crj,o = Im/s, i.e. high values are penalized. For a discussion of 
priors, see iFord and Gregorvl (j2007| ). We adopt a likelihood that is the product of two terms corresponding to the 
radial velocity and transit observations. The likelihood for radial velocity terms assumes that each radial velocity 
observation ivi) is independent and normally distributed about the true radial velocity with a variance of a1 + (t|, 
where is the published measurement uncertainty. 

Instead of modeling each photometric observation, we account for the transit observations by including a likelihood 
term that is the product of three Bayesian penalties based on the orbital period, transit duration and ingress time of 
HAT-P-13b, assuming Gaussian distributions for each measurement with standard deviations taken from t he pub lished 
uncertainties as derived by B09. We use MCMC to calculate a sample from the posterior distribution (jFordl [20061. 
We calculate multiple Markov chains, each with ~ 2 x 10^ states. We discard the first half of the chains and calculate 
Gelman-Rubin test statistics for each model parameter and several ancillary variables. We find no indications of 
non-convergence. Thus, we randomly choose a subsample (10, 000 samples, large enough to give a statistically valid 
and accurate outcome, but small enough to be computationally tractable) from the posterior distribution for further 
investigation. 

2.2. Transit Timing Variations 

To investi g ate th e T TV signature in th e HAT-P-13 system, we employ the same basic method as that used in 
iPavne et al.l ()2010f) and iVeras et al.l ()2010f ). We consider a fiducial system consisting of two planets: a transiting 
hot- Jupiter planet and an outer, non-transiting planet which perturbs the transit times of the inner planet. Given 
an initial specification of the planetary masses and orbital elements, we evolve each system forward in time for ~ 3.5 
years, corresponding to several hundred transits by the inner planet in a system such as HAT-P-13 (assuming that 
the inner plan et remains tran siting throughout the in tegration). This 3.5 year integration time was also used in the 
studies of Pav ne et al.l ()2010D andj^ras et al] (I2O1O0 whose investigations attempted to illuminate potential Kepler 
missions observations, where the Kepler mission is expected to run for at least 3.5 years. To facilitate any comparison 
with the method and results of these papers, we chose to also maintain this 3.5 year simulation timescale. 

The n-body integrations are performed using a conservative Bulirsch-Stoer integrator, derived from that of Mercury 
(Chambers 1999). We use a barycentric coordinate system and limit the time steps to no more than 0.04 times the 
orbital period. After each time step, we test whether the star-planet separation projected onto the sky (A) passed 
through a local minimum and the planet in question is closer to the observer than the star. Each time these conditions 
are met, we find the nearby time that minimizes A via Newton- Raphson iteration and increment an index i. If the 
minimum A is less than the stellar radius, then we record the mid-time of the transit, ti. In calculating the observed 
transit time one needs to account for the light travel time, S titt(i) — — (rp(ii) • ^ios)/c, where Vpiti) is the barycentric 
vector of the planet at time ti, rjos is a unit vector pointing to the observer, and c is the speed of light. The observable 
transit time variations are calculated as S t{i) = ti -\- S tiu{i) — i P — to, where the constants P and t^ are determined 
by linear least squares minimization of '^^{5 t{i))'^. We neglect any motion of the stellar center between the time of 
light emission and the time of transit. 

It should be noted that the TTV investigations in this paper split into two main strands: (i) an investigation in 
N3.2I of the TTVs that would arise in the systems that come out of our RV MCMC analysis, and (ii) the structured 
investigation of inclination effects in TTVs for HAT-P-13 carried out in ijH The introduction to the respective results 
sections provides more detail on the precise manner in which each TTV investigation was conducted. 
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3. MCMC ANALYSIS OF THE RADIAL VELOCITY OBSERVATIONS 

3.1. Jitter 

In the B09 diseovery paper, upon performing an RV fit (in conjunction with an analysis of the observed transit timing 
constraints) they found that a jitter of 3.0 ms~^ was required in order to give a reduced figure of 1.0. Similarly 
WIO found a jitter of 3.4 ms^^ was required in order to give a reduced figure of 1.0. In the work of lWrightl §005), 
it was found that from a sample of ^ 30 stars similar to HAT-P-13 (mass, Ah = 1.22Mq) that the distribution of 
jitters was such that the 20th, 50th and 80th percentiles were found to be 2.6, 4.0 and 6.2tos~^ respectively, with 
an upper bound of ~ 15ms~^. It is thus plausible that the actual jitter value for HAT-P-13 is substantially different 
from the value of ~ 3.0 ms~^ used in B09 and WIO. As discussed in the introduction, we feel that it is important 
that the jitter be placed on an equal footing and analyzed in the same manner as all of the other parameters that one 
can constrain using the RV data. As such, we include the jitter as a prior in our MCMC analysis (see Section [^?T] for 
details). We illustrate in Fig [T] the sensitivity of a number of the key parameters to the value of the jitter that is used 
in the MCMC analysis. 

To illustrate the general effect of including jitter as an integral part the MCMC analysis, we look at two different 
data cuts for IIAT-P-13: The first set (SI) uses only the RV data that was included in the discovery paper of B09, 
while the second set (S2) uses the larger data set used in the later paper of WIO. We have analyzed these data sets in 
such a manner as to afford direct comparison with the original papers: we analyze SI assuming a two-planet solution 
(as was done in B09), while we analyze S2 by allowing for the possibility of an additional "slope" in the RV data 
(see Eqn 1 of WIO) Additionally, for S2 we also perform a full 3-planet analysis. We plot some of the results of these 
MCMC analyses in FiglH 

In general, it can be observed from the various plots in Fig [1] that the MCMC routine investigates a wide range 
of jitter levels, and that (as one might expect), the larger the level of jitter it tries, the wider the range of planetary 
parameters it can accommodate. However, we see that this process is not completely unlimited, as there are few points 
in any of the plots which have jitters above 30 ms""'^, and the majority of them are significantly below this level. 

In more detail, we see that from the left-hand column of Fig [T] (in which the SI results are plotted), that if we 
restrict ourselves to examining jitter ranges 2.0 < aj < 3.0tos~^ (see the green crosses (gray in the print version) 
in Fig [1] which label regions with aj — 3.0 ms~^ ± 50%), the eccentricity of the inner planet, emner, is constrained 
to be O < einner < 0.07 (best fit value of ei„„er = 0-017lo;Qjg, where the uncertainty figures cover 68% of the 
data, i.e. they equate to 1-sigma error bars), whilst the difference in the alignment of the longitude of pericenters, 
77 = Wouter — Thinner, Can take ou essentially the full range of available values: —180° < t] < 4-180° (best fit value of 
r/ = 15.9^°g3i^). If we allow a higher range of possible jitter, then the available parameter space expands, such that if 
we take all of the MCMC results with < aj < lOms^^, the eccentricity of the inner planet could take any value in 
the range < einner < 0.15 (best fit value of einner — 0-021 19012)' while the eccentricity, eouter, and period. Pouter, 
of the outer planet also expand to encompass ranges of 0.6 < eouter < 0.95 (best fit value of eouter — O-^^^om) ^^-^ 
410 < Pouter < 450 days (best fit value of Pouter = 430.2^4 g) respectively. Interestingly, we note that at such high 
values of aj, the alignment of the longitude of pericenters starts to favor orthogonal values, rj « ±90°. 

The middle column in Fig[l]details the results of the 2-planet -I- slope analysis of S2, i.e. analogous to the analysis of 
WIO. We see that there are now very few systems in the analysis which have jitter levels below 5 m s^^. If we restrict 
our attention to the very few systems which have aj — 3.0m ± 50%), the eccentricity of the inner planet, einner, 
is constrained to be < einner < 0.09 (best fit value of einner = O.OSSIIq'^H, a value very much larger than that found 
in SI (or indeed, in the analysis of WIO). If we allow a higher range of possible jitter, then the available parameter 
space expands, such that if we take all of the MCMC results with < CTj < 10 ms~^, the eccentricity of the inner 
planet could take any value in the range < einner < 0.2 (best fit value of einner = 0.06j;Q 026)' 'while the eccentricity, 
Bouter, and period. Pouter, of the outer planet also expand to encompass ranges of 0.45 < eouter < 0.9 (best fit value 
of 

^outer 

— 0.64to;Jg) and 430 < Pouter < 490 days (best fit value of Pouter — 450.617^) respectively. 

Finally, in the third column of FigdJ we display the results from our 3-planet analysis of S2. The results pertaining to 
Sinner, ^outer. Pouter and 7] as displayed are qualitatively very similar to those of the 2-planet -I- slope analysis, as find 
that for CTj = 3.0 ms^^ ± 50% the best fit values are einner = 0.087lQ;o|g, eouter — ^-^^^^'m ^^"^ Pouter — 455.0ti3.7. 
The fitted parameters for the third planet in this analysis are extremely ill-constrained. Given the qualitative similarity 
in the two sets of results for the S2 system for the inner planets and the very poor constraints that can be placed on 
any third planet, we chose to concentrate the rest of the analysis on the 2-planet + slope analysis to ensure ease of 
comparison with the published results of WIO. 

The goodness-of-fit parameters plotted in the final row are an effective-chi-square measure, used here to illustrate 
the rather large range of fit "qualities" that the different MCMC runs result in. We can see that in the 2-planet 
analysis of B09 in the left-hand column, it is particularly noticeable that the low jitter systems tend to have a worse 
fit (higher effective chi-square measure) than do the systems with higher jitter. The analyses of the WIO data show a 
wide range of effective chi-square measures right across the jitter range sampled by the MCMC routines. 

We provide two tables at the end of the subsequent section (Table [Hand Table [2|) in which we list some key statistics 
for the fitted values of the various system parameters arising from our respective analyses of SI and S2. We give the 
median value for each parameter, and then also give as uncertainty figures the spread in parameter values which are 
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Fig. 1. — Sensitivity of selected system parameter determinations to assumptions regarding jitter. The top row has Souter vs. Jitter, the 
second row has Pouter vs. Jitter, the third row has rj = zuouter — dinner vs. Jitter, and the bottom row has Goodness-of-fit vs. Jitter. 
In the left-hand column we present results using only the subset of data known at the time of publication by B09, assuming a 2-planet 
fit. In the center and right-hand columns we present results obtained using the full data set of WIO, with the central column assuming 
2-planets -|- a linear trend, while the right-hand column assumes a 3-planet fit. Systems with aj < 2.0m are plotted using a red filled 
circle (gray in the print version), those with 2.0 < aj < 4.5 ms~^ are plotted using a green cross (gray in the print version), and those 
with cTj > 4.5 ms~^ are plotted using black -|- symbols. If one regards the jitter as being unknown and allow a range of values to be 
sampled by the MCMC routine, we can see that in this sample of plots for various system parameters, there are a much larger range of 



parameters which can plausibly fit the data than if one assumes CTj 
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required to cover 68% of the data in the sample, i.e. corresponding approximately to the 1-sigma deviation figure as 

quoted in B09. 

We discuss further some of the implications of these jitter figures in i )3.3l 

3.2. Combining RV and TTV Data 

In their discovery paper, B09 also perform an analysis of the transit timings to look for any evidence of TTVs. They 
find that the TTVs in the system are restricted to be < 100 seconds. We wish to add-in this constraint directly to the 
MCMC analysis performed in section [3l To do this, we take a 10,000-strong subsample of the MCMC systems (see 
section [5] for further discussion of methodology) and use these systems as the basis for a TTV investigation. 

We take the fitted parameters for each of the subsample systems, and use these to set up the planetary masses and 
orbits in an n-body simulation. Given that the relative inclination of the planets is unconstrained in the RV analysis, 
we performed the simulations assuming coplanarity. This simulation is then run in the manner described in section [2?2] 
and a TTV analysis performed on the inner planet. This allows us to directly investigate whether any particular range 
of the MCMC subsample can be excluded by considering the observational TTV constraints. Extremely high TTV 
signals generally arise due to close approaches between objects in the simulation, giving a strong indiction that such 
systems are unstable (see Veras et al 2010 and Payne et al. 2010 for further discussion of such high TTV signals, and 
the likely Hill and/or Lagrange stability.). However, irrespective of whether any such systems are absolutely stable or 
not, we can use the fact that they have TTV amplitudes > 100s to remove them from the analysis. 

We show a sample of the results from the TTV analysis in Fig[2j On the left-hand side we plot results for a 2-planet 
analysis of the SI data set, while on the right-hand side we plot results for the 2-planet -I- slope analysis of the S2 data 
set. The majority of the observable parameters do not show any obvious correlation with the RMS TTV amplitude 
(we provide some examples of such plots in an appendix). However, the overall RV offset, as well as the RV Amplitude, 
Kouter , longitude of pericenter, Wouter, and the eccentricity, Couter, of the outer planet all have such a correlation. We 
could thus hope to use this upper TTV amplitude cut-off of 100s to limit the allowed range of these variables. Whilst 
this is certainly possible on a gross scale, we note that when we color-code the results according to the assumed jitter 
(as done in Fig[2|), the majority of the low jitter systems tend to fall into the low-TTV area of parameter space. This 
means that using the TTV amplitude constraints is probably only of some marginal help in restricting the range of 
plausible system parameters to be considered. 

In Fig [3l we go on to plot a selection of the key orbital parameters, using colored symbols to show the regions of 
parameter space which are favored or ruled out by the combined jitter and TTV amplitude analyses discussed above. 
As in Fig [21 on the left-hand side we plot results for 2-planet analysis of the SI data set, while on the right-hand side 
we plot results for the 2-planet -I- slope analysis of the S2 data set. These suggest that for the inner planet, even for 
a fairly tightly constrained jitter of 2.0 — 4.5m s~^, there can be substantial variation in the eccentricity, einner, such 
that it can take any value < Cinner < 0.07; 

For the outer planet, much greater variations are possible: (i) Perhaps the strongest result we find from using an 
observed constraint on the TTV amplitude of < 100 seconds, is in constraining the eccentricity of the outer planet 
to be Couter < 0.85, irrespective of assumed jitter values, giving an approximate overall range of 0.5 < Couter < 0.85 
(best fit value of Couter = O-OlT^QQQg, assuming jitter is in the range 2.0 to 4.5 ms~^). This is clearly comprehensible 
in the light of the above argument concerning stability: extremely high eccentricity values for the outer planet can 
easily lead to close approaches and/or orbit-crossing. These results stand for the analyses of both SI and S2. (ii) 
While the longitude of pericenter of the outer planet (wouter) is strongly constrained to be ~ 180 (best fit value 
of Wouter = 177.0lQg for SI, TUouter = 171.2j^y'3 for S2), the longitude of pericentre of the inner planet can take 
almost any value, meaning that there is no strong preference to indicate that the planetary pericentres are aligned 
(-180° <rj< +180°, best fit value for SI of 77 = Id.OtH l, '^hile for S2 the best fit value has rj = 3.21^^;^). (in) In the 
SI analysis, the outer planet can have periods 420 < Pouter < 440 days (best fit value of Pouter = 429.7l3 g), while in 

the S2 analysis, the period is pushed to higher values 430 < Pouter < 490 days (best fit value of Pouter = ^^b.Otif^)- 
Note that Table [T] also contains a section in which the figures are analyzed for systems with TTVs less than 100 
seconds only. As an example of the data contained therein, it can be seen that the lower allowed range on the outer 
planet's eccentricity from the TTV analysis manifests itself in the table as a reduction in both the upper uncertainty 
and the median value of the fitted outer eccentricity. 

3.3. Summary of MCMC Results 

Our re-analysis of the radial velocity data and subsequent introduction of a coupled TTV analysis leads to the 
following conclusions: 

Comparing our results with those of the HAT-P-13 discovery paper B09 (where aj — 3.0m was used), we find 
that a rather larger range of fits to the RV data are possible. In particular, we find (for the SI data set used in B09) 
that even a rather modest levels of jitter (2.0 < aj ^ 4.5 ms~^) wiU allow Cinner = 0-017to.oo9J Pouter = 0.73j;Q |!)g, 
Pouter — 429.7^3'7 days, and relative pericenter alignment, rj = 15.9'^^^'^. 

Similarly, when we look at the S2 data set as used in WIO, a level of jitter 2.0 < aj ~ 4.5 rns^^ will allow 

einner = 0.038^0^11, Router = 0.65lo;Jg, Pouter = 449. 4t^;^ days, and relative pericenter alignment, rj = -12.4t|3;9. 

Table [1] contains a summary of all the parameters and the best fits resulting from our analysis of SI, while Table [2] 
contains a summary of all the parameters and the best fits resulting from our 2-planet -|- slope analysis of S2. We note 
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Fig. 2. — Understanding the constraints given by combining the TTVs and the Jitter: I. TTV Amplitude vs. System Parameters. In 
the various plots above we show the RMS TTV amplitude as a function of various system parameters (planetary eccentricity, etc). As is 
Fig. 1, systems with uj < 2.0 ms~^ are plotted using a red circle (gray in the print version), those with 2.0 < < 4:.5ms~^ are plotted 
using a green X (gray in the print version), and those with aj > 4.5m are plotted using black + symbols. In the left-hand column we 
present results using only the subset of data known at the time of publication by B09, assuming a 2-planet fit. In the right-hand column 
we present results obtained using the full data set of WIO, assuming 2-planets -f a linear trend. We can see that imposing a cut such that 
the TTV amplitude is < 100s helps to eliminate a sizeable proportion of the overall solutions, but note that the majority of the systems 
with jitter < 3.5 ms~^ tend to naturally fall into areas which have TTV amplitudes < 100s. The quantities plotted above are the only 
variables which have a potentially significant correlation with the RMS TTV amplitude. We note that comparing the WIO data to the 
B09 data, once again the TTVs can be used to (approximately) constrain the eccentricity of the outer planet, but now is even worse at 
constraining the other parameters due to the much larger scatter. 
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Fig. 3. — Understanding the constraints given by combining the TTVs and the Jitter: II. We plot a variety of the possible system 
configurations arising from the MCMC analysis, and then constrain the possibilities by indicating (i) Using blue squares (gray in the print 
version) the systems which have TTV amplitudes bigger than 100s, and then (ii) In other colors indicate the range of jitter assumed for 
the systems: systems with < 2.0ms~^ are plotted using a red circle (gray in the print version), those with 2.0 < crj < 4.5ms~^ are 
plotted using a green cross (gray in the print version), and those with aj > 4.5m are plotted using black addition symbols. As in Fig. 
1, in the left-hand column we present results using only the subset of data known at the time of publication by B09, assuming a 2-planet 
fit. In the right-hand column we present results obtained using the full data set of WIG, assuming 2-planets + a linear trend. We find that 
the eccentricity of the outer planet is strongly constrained to be < 0.85, but that other parameters are largely unaffected / unconstrained. 



that within the range of our uncertainties, all of our results are consistent with those of the discovery paper B09, but 
we do in general find a much broader spread of allowed values for all parameters. Similarly for the comparison of the 
S2 analysis with the results of WIO, our results are consistent with those of WIO (except for the inner eccentricity) but 
the error bars on our fitted parameters are much broader. To "make the eccentricities consistent" , one has to appeal 
to even lower jitter values 2.0 < Cj ~ 3.5 ms~^, by which point the lower 1-sigma limits from our analysis start to 
overlap with the upper limits from the WIO analysis. 

We should point out that while our prior on the jitter (§2.1) does penalize high jitter values, it may be that an even 
more punitive prior is warranted. This is due to the observed fact that stellar activity (and therefore stellar jitter) 
declines with age (,Isaacson and Fischer .2010.) . As HAT-P-13 is old (~ 5Gyr, it is possible that it has a jitter level 
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TABLE 1 

W e provide a table summarizing the results of our MCMC analysis for the HAT-P-13 system using just the data from the work 
of \Bakos et al.\ \2009i)- The different columns give the figures which result from different jitter assumptions, with the top 

HALF OF THE TABLE SHOWING RESULTS FOR all DATA, WHEREAS THE BOTTOM HALF SHOWS THE DATA FOR JUST THOSE SYSTEMS FOUND TO 
HAVE TTV < 100 SECONDS. THE FIGURES GIVEN ARE THE median VALUES FROM THE 10,000 SELECTED MCMC REALIZATION, WITH THE ± 
FIGURES GIVING THE SPREAD REQUIRED IN ORDER TO ENCOMPASS 68% OF THE DATA, I.E. GIVING A FIGURE WHICH APPROXIMATES THE 
1-SIGMA RESULT TABULATED IN THE DISCOVERY PAPER Of IBAKOS ET AL.I II2009I) . 

TTV Parameter a^^ < 3.0 < 4.5 a^^ < 10.0 All a'^ 2.5 < a'^ < 3.5 2.0 < a'^ < 4.5" 



ALL 


dinner (days) 


2 qi626+°"°''°2 


2 91626+0 00002 

Z.310Z0_Q gggg2 


2 91626+0 00002 

Z.310Z0_g gggg2 


2 91626+0 00002 

Z.31UZO_g gggg2 


2 91626+0 00002 
^■^^"^"-0.00001 


2 91626+0-00002 
^•^l"^"-0. 00002 


ALL 


dinner 


Q Q1 7+0.009 

U.UH_Q^QQg 


017+0-012 

U-Ui/_Q QQg 


0.021^0.023 


0.023^0.032 


016+0.012 

U.U1D_Q Qgg 


017+0-013 

U-U1(_Q ggg 


ALL 


inner 


157.1_37 g 




164.9l^^j 






1C1 n+86.1 

161-9_go.i 


ALL 


P, inner 


4780. 5lg;^ 


4780.61q5 


4780. 81q 7 


4780.91q g 


4780. 4 


4780. etgg 


ALL 


Pouter (days) 


429.4+|-^ 


429.7+^-* 




430.3+^-^ 


429.9+3 g 


429.7+^* 


ALL 


^outer 


'J- '^-0.04 


73+0.12 

"•'■'-0.06 


75+0-12 
'^•'^-0.08 


75+0.12 
'^•'^-0.08 


72+0.12 


7S+0-12 


ALL 


outer 


177 1+0-^ 


1 77 2+l'0 


177.4ti;2 


177.4^1-4 


177 1+1-0 

1' '-1-0.9 


177.2ti;o 


ALL 


P, outer 


4890. l+g® 


4890.2t0-7 


4890. ll?o 


4890.lti-2 


4890. lto7 


4890.2+0-^ 


ALL 


V 


20.6l?«i 




ir, + 64.8 

-^^■'^-95.8 


ir, 9 + 67. i 

.^^■^-98.5 


-17 7+44.9 
l'-'-81.6 


15.91^3:? 



ALL 


Tt, outer, Bakos"^ 


4872.4I4 


4874.0j;g;5 


4875.9t7;g 


4876.lt^;J 


4873.0t4;3 


4874.lt^;2 


ALL 


Tt, outer, 2010'' 


5302.6tg:g 


5304.6l7:g 


5306.7lg''g2 


5307.0tio4 


5303.9^7:3 


5304.6l7:g 


ALL 


TT,outer, 2011^ 


5732.1+9^2° 


5734.3^^2-1 


5736.7^12 8 


5737.1^13 1 


5733.8^^2-7 


5734.3lig-2 


ALL 


Tt, outer, 2012^ 


6161.3^1^-0 


6164.0^15-] 


6166.6t2i.;2 


6167.lt?|-;* 


6164.1^15-6 


6164.lti^-| 


< 100s 


Pinner (days) 


2 91626+0-00002 
^•^l"^"-0. 00001 


2 91626+0-00002 


2 91626+0-00002 
z.yiuzu_Q 00002 


2 91626+0-00002 
z.aiozo_Q 00002 


2 91626+0-00002 
z.yiozo_o 00001 


2.91626tO:S2 


< 100s 


dinner 


Q (,,7 + 0.009 


017+0-013 
U-Ul/_o 009 


021+0-023 
U.UZ1_Q 012 


0.023to«f3 


016+0-012 

U.U1D_Q 009 


017+0-013 
U-Ul/_o 009 


< 100s 


inner 


155.8_3g;4 


161.21^:? 


1CK Q + 96.0 
165.3_62.7 


165.5t^^:« 


156.8l^5:g 


iDi.U_4g g 


< 100s 


P, inner 


4780.5^15:3 


4780.6lo 4 


4780. 8lo;? 


4780.9lo:g 


4780. 5loi 


4780.6lo:5 


< 100s 


Pouter (days) 


429.4+i^-i 


429.7111 


430.2^4^ 


430.3^5 2 


430.0+3 9 


429.7+^1 


< 100s 


^outer 


70+0-07 
U.'U_g 134 


71+0.07 
'J-'i-o.os 


p, 79+0.08 
'^-0.06 


p. 79+0.08 
'^-0.06 


71+0-07 
'^•'1-0.04 


72+0-07 
'J- '^-0.05 


< 100s 


'^outer 


176.9tg;^ 


177.0l°i 


177.0t?;? 


177.0^^!^ 


176.9tg? 


177.0t0-i 


< 100s 


rpa. 

PfOuter 


4890. 1+07 


4890. 1+og 


4890. ll?o 


4890. Olij 


4890. 1+0 7 


4890. llog 


< 100s 


V 


20.6l?^:| 




7+64.2 
11- '-94.6 


ii-'^-oe.g 


20.5tg:| 


16 0+50.2 

iD.U_g2 g 



< 100s Tt, outer, Bakos'^ 4871. 7t3 g 4872.514 5 4873. itg g 4873.ltg i 4872.oj;3 g 4872. et^ j 

< 100s TT,o„ter,2010^ 5301.5lg:| 5302. 5lg 4 5303. ell;® 5303.7lg;5 5301.9lg:g 5302.5lg:5 

< 100s Tt, outer, 2011'^ 5730.5lg''o5 5732.2lgi40 5733.5tii;4 5733.7^12.4 5732.01922" 5732.2lgi4^ 

< 100s TT,o„ter,2012=' 6159.9tiof 6161.7ti2g 6163.5^20.1 6163.8^21-8 6162.oti3o 6161.7^125 
We strove to make our data comparable to the discovery paper where at all possible: as such, it should be noted that the parameter fits 

given for P, e, ro and Tp^ri are based on a coplanar MCMC analysis, while the transit time calculations for the outer planet are made 
using the assumption that the outer planet is inclined at 90° to the plane of the sky. 
^ All transit and pericenter passage times are in Julian days - 2,450,000 
1" All jitter constraints are in ms~l 
° All pericenter alignments are in degrees 



that is drawn from a population with lower values than that quoted from the general population in I Wright (|2005h . As 
such, it is possible that a more limited range of jitters should be allowed. However, we note that (i) our jitter allows 
for contributions over and above pure "stellar jitter" , and (ii) we list results for a number of different jitter ranges in 
Tables [T] & [2] from which we see that even if one only wishes to accept a limited ranges for the jitter values resulting 
from the MCMC analysis (e.g. 2.0 < aj ~ 4.5 ms~^), the associated errors on the fitted planetary parameters are still 
much broader than those presented in the results of either B09 or WIO. 

Given that any "jitter" values arising from our analysis will intrinsically be "system" jitter (the combined contribution 
of jitter from stellar sources, instrumental noise, undetected planets, etc) rather than pure stellar jitter, we are hesitant 
to produce detailed figures regarding any "best-fit" jitter parameters arising from our analysis. However, for the same 
of completion, we note that the analyses corresponding to the three columns of Fig 1 (SI, S2 assuming 2 planets + 
slope, and S3 assuming 3 planets) find overall system jitter levels of Uj = , aj = and aj = X^^ respectively. 

Subsequent to the pure MCMC analysis, if we then introduce an additional TTV analysis step (based upon the 
output of the MCMC routine), we can somewhat constrain the eccentricity of the outer planet, finding that now 
0.6 < Couter < 0.85. However, the remainder of the other elements are essentially unconstrained by this additional 
analysis. Tables [T] & [5] also contain a summary of all the parameters and the best fits resulting from our MCMC + 
TTV analysis. 

We thus conclude this section of our investigation by suggesting that previous methodologies of excluding jitter 
from an MCMC analysis significantly overestimates the precision with which other system parameters (e.g. planetary 
orbital elements) can be estimated from RV data. We thus recommend that jitter be included as a model parameter 
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TABLE 2 

We provide a table summarizing the results of our MCMC analysis for the HAT-P-13 system using the full set of data 
available from WIO. The different columns give the figures which result from different jitter levels, with the top half of 

THE table showing RESULTS FOR all DATA, WHEREAS THE BOTTOM HALF SHOWS THE DATA FOR JUST THOSE SYSTEMS FOUND TO HAVE TTV 

< 100 SECONDS. The figures given are the median values from the 10,000 selected MCMC realization, with the ± figures 

GIVING THE SPREAD REQUIRED IN ORDER TO ENCOMPASS 68% OF THE DATA, I.E. GIVING A FIGURE WHICH APPROXIMATES TO A 1-SIGMA 

RESULT. 
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0.03810;°- 
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-0.4 
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-0.06 
5 + 1.9 



449, 
0.65 
175.3 
23.3 



+45.3 
82.9 



+0.00004 
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We strove to make our data comparable to the discovery paper where at all possible: as such, it should be noted that the parameter fits 
given for P, e, zu and Pperi are based on a coplanar MCMC analysis, while the transit time calculations for the outer planet are made 
using the assumption that the outer planet is inclined at 90° to the plane of the sky. 
^ All transit and pericenter passage times are in Julian days - 2,450,000 
All jitter constraints are in ms~^ 
° All pericenter alignments are in degrees 

in future MCMC analyses of RV data, allowing a true estimation of system parameters and their uncertainties to be 
determined. 



4. FURTHER INVESTIGATION OF TRANSIT TIMING VARIATIONS 

In the coplanar MCMC -f TTV analysis of section [31 it was suggested that the RMS TTV amplitude could be of 
so me diagno s tic po wer in constraining cert ain system parameters (E.g. the eccentricity of the outer planet, Couter)- 
In iNesvornvl (|2009f ) and lPavne et al.l ()2010[ ). it has been found that planets which are significantly inclined relative to 
one another can give rather different signal profiles and amplitudes as compared to the coplanar case. This inclination 
dependence, combined with the differing planetary masses that would be required in an inclined system to satisfy the 
RV constraints, suggests that a more detailed investigation of the inclination dependence of the TTV signal in the 
HAT-P-IS case is warranted. 

Our ultimate aim is to find whether it might then be possible to use the TTV characteristics to break the degeneracy 
of the RV M outer sin i outer = 15.2Af / relation. 

4.1. Examples of Transit Timing Variation Signals for HAT-P-13 

We commence this more detailed investigation of the inclination dependence of the TTV signal variation in the 
HAT-P-13 system by plotting in FigUsome sample results which illustrate the effects of varying the orbital inclination 

of the outer planet. Note that we set Oinner = 0.043 At/, Cmner = 0.02, rUinner = O.SSM,/, Uouter = 1.188, eouter = 

0.691, &z ruouter siniouter = 15.2Mj in line with the standard values from B09. N.B., as we vary the inclination of the 
outer planet, iguter, we also simultaneously vary the mass of the outer planet, niouter- 
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Fig. 4. — Variation of the relative planetary inclination in the HAT-P-13 system. The system parameters are as observed, such that 
o-innur = 0.043 Ai7, Cinner = 0.02, minn&r = 0.85Mj, Uouter = 1.188, Couter = 0.691, niouter siniouter = 15.2Mj. While the pattern of 
behaviour is non-trivial, we can make certain statements: (i) For the most likely prograde cases i ~ 0° and i ~ 45°, there is a significant 
difference in the predicted TTV profile, with much sharper features being expected for the high inclination cases; (ii) The TTV amplitude 
for these cases is ~ 20 seconds, well within the current observational constraints; (iii) Systems with planets close to perpendicular are likely 
to give highly distinctive, high amplitude TTV signals; (iv) The retrograde systems can give rather similar amplitudes and overall profiles 
to those of systems in prograde orbits (extra care may be needed to analyses and distinguish these cases). It should be noted that the 
origin of the time axis on this plot is somewhat arbitrary, given that it is a simulation, but the translation to "real-life" can easily be done 
by noting that the deepest troughs in the TTV plots occur as the outer planet approaches pericenter, e.g. April 16th, 2010. 
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Fig. 5. — TTVs as a function of time for various values of the outer planet eccentricity, eouter, in the prograde co-planar HAT 13 b 
system. We find from this particular example that increasing eouter from 0.5 to 0.8 causes only a minor increase in the TTV amplitude, 
but that going to higher eccentricities (eouter ^ 0.85) leads to the significant amplitude variations seen in Fig [21 potentially associated 
with unstable planetary orbits. As in Fig (4] it should be noted that the deepest troughs in the TTV plots occur at the time of pericenter 
passage of the outer planet. 

Given the work of lMardlin^ ([2010') discussed in gU we concentrate our analysis on a few key relative inclinations; 
For the prograde cases we concentrate on irei ~ 0° and irei ^ 45°, as these are thought to be the most likely states at 
late times. We also examine the symmetric retrograde cases, irei ^ 180° and irei ~ 135°. Finally, for completeness we 
also look at a couple of example plots from the "forbidden" region of parameter space 55° < irei < 125° from which 
it is thought the planets are likely to be excluded due to the combined effects of Kozai forcing and tidal dissipation. 

We can immediately see from the irei = 0° and irei = 45° plots in Fig HI that (a) the expected TTV amplitudes 
are well below the current observational limits, so either inclination is plausible in that sense, but (b) perhaps more 
importantly, the profiles have rather different shapes, with the irei = 45° plot exhibiting prominent spikes in TTVs 
over a relatively short period of time, in addition to having an overall reduced amplitude. We discuss further in section 
[S] a plausible observational strategy to extract this information. 

We can also see that the completely retrograde case, irei = 180°, is very similar in both shape and amplitude to the 
irei — 0° prograde case, suggesting that they would be rather difficult to distinguish observationally. 

For the plots in the "forbidden" 55° < irei < 125° region, it seems that the systems in which the planets are close to 
perpendicular will have extremely high TTV amplitudes (due, no doubt, to the significantly increased iriouter required 
in order to satisfy the mouter sin? constraint). These highly distinctive profiles mean that if for some reason a planet 
were able to occupy this region of parameter space, it would be readily identifiable (and indeed, certain particularly 
high inclinations may already be ruled out by the TTV constraint - see section 14.21 below for further discussion) . 

Finally, we note in passing that in section |3] it was found that eccentricities for the outer planet were constrained to 
lie below eouter ^ 0.85. We show in Fig[5]some details of the TTV plots which result from increasing the eccentricity of 
the outer planet. Keeping the semimajor axes and masses of the planets fixed at the best-fit values from the discovery 
paper (and with irei = 0°), we increase eouter from 0.5 to 0.85. We find that for values of the eccentricity in the 
range 0.5 < eouter < 0.8, increases in eccentricity result in only marginal increases in the TTV amplitude. However, 
at very high eccentricities {eouter ~ 0.85), the systems start to exhibit "kicks / step-changes" in their behavior as 
the outer planet starts to significantly perturb the inner orbit (as the outer planet passes through pericenter). We 
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note that whilst such systems may be stable against collisions, the relatively large perturbations experienced during 
close pericentre en count ers can still l ead to significant orbital evolution and hence large TTV amplitude variation (see 
iPavne et all (|20Tol) and lVeras et"all ^OW) for more detailed discussions). 

4.2. Transit Timing Variation Contour Maps 

To gain a further insight into the range of TTV signatures possible for HAT-P-13, we keep the mass and orbit of the 
inner planet fixed and vary the orbit of the more massive outer planet. For each set of configurations - i.e. semi-major 
axis, a, eccentricity, e, & inclination, i of the non-transiting planet - we typically conduct 5 simulations randomising 
the angular orbital elements (argument of pericentre, w, longitude of ascending node, i7, & initial mean anomaly M). 
We assume that the known transiting planet has the mass, semi-major axis and eccentricity of HAT-P-13b {Minner — 
0.85 Mj, ainner = 0-043 AU and Cmner = 0.02) and constrain the outer planet to have M outer sin i outer — 15.2 Mj). 

In Fig |6] we calculate the RMS of 5t{i) for each configuration and present contour plots of the median RMS TTV 
amplitude in the ( '^outsr ^ Couter) plane for multiple inclinations. Note that these plots use simulations conducted for 
3.5 year simulations, but similar simulations conducted over 10 year periods are essentially identical. Moreover, given 
that the plots showing the maximum RMS TTV amplitude and those showing the median RMS TTV amplitude have 
qualitatively and quantitatively similar values across the parameter space, we use the median, 3.5 year plots for all 
figures below. 

Looking at the irei = 0° prograde case in Fig |6l we find that the values of the outer planet's semi-major axis and 
eccentricity as given in the B09 or WIO papers(ao„ter ~ 1.2 AU & Couter ~ 0.6 — 0.7) would be expected to give an 
RMS TTV amplitude of ~ 10s. Furthermore, we find that there is a rather large range of parameter space which gives 
similarly low TTVs, meaning that the current TTV observations cannot help to constrain the outer orbit, other than 
saying that Couter < 0.8 {oouter is essentially unconstrained by the TTVs in the range 1 < aouter < 1.5 AU which was 
examined) . 

As the relative inclination of the orbits is increased, we find that there is a slight decrease in TTV signal amplitude 
in the Couter ^0.7 region as i — 45 °. 

As the relative inclination of the orbits increases further towards irei = 90 °, there is an increase in RMS TTV signal 
amplitude across the entirety of the parameter space, such that at irei = 85 ° almost all regions have amplitudes > 30 
seconds. For i = 90°, the mass of the outer planet would become formally infinite, so we cannot investigate this case, 
but we do note that in the irei = 89 ° case, almost the entirety of parameter space has TTVs higher than the observed 
limits. 

As the orbits become retrograde, we find that the TTV amplitudes decrease again, with the results at irei — 135 ° 
(irei = 180°) being very similar to those seen at irei =45° {irei = 0°). I.e. there is some kind of approximate 
symmetry about Zrei = 90°. 

The clear difference in TTV amplitude as i 90 ° is a rather different behaviour to that observed by iPayne et al.l 
()2010( ) in their general investigation of inclinations effects on TTVs in hot- Jupiter + hot-Earth systems. There they 
found that TTV amplitudes were very generally highest at irei = ° , with an approximately monotonic decline as 
i — >■ 180°, with no special behaviour se en at irei =90°- H owever, it is difficult to make a precise comparison between 
their results and our results here, as in IPavne et al.l ()2010[ ) the planetary masses were kept constant as the inclination 
was increased, whereas in this investigation we are varying the outer mass as Mouter = 15.2 Mj/ siniouter- 

For completeness, we include in Fig [7] a plot in the {irei, Pouter) plane, again plotting the median RMS TTV values 
from 3.5 year simulations. In these simulations we fix the semi-major axis of the outer planet and then vary the 
inclination and eccentricity of the outer planet. We reproduce here the plot for aouter / o-inner = 27.2, i.e. an outer 
planet of Pouter = 429 days. As seen from Fig [71 the upper allowed eccentricity decreases at larger semi-major 
axes. There is a critical range of inclinations above which the TTV amplitude becomes too great to fit observational 
constraints: For the current 100 second limits, the excluded range is very narrow ~ 90° ± 2°. Clearly this is of little 
benefit, as such regions are already excluded in any practical sense, as such a highly inclined object would be so massive 
as to be a star and hence likely be readily visible, but if future observations constrain the TTV amplitude to be less 
than 10s (for example), then the inclination restrictions become much more severe, and all inclinations ~ 90° ± 20° 
could be ruled out. Moreover, we see that as the inclination increases from 0° — 45°, the ability to constrain the 
eccentricity using TTV limits varies significantly: E.g. the yellow region mapping out the 10 — 30s TTV contour 
narrows rapidly. 

We also made similar plots at different semimajor axes, but we find that there is little difference between the results 
at different semimajor axes, as might be anticipated from Fig|6l 

5. FUTURE TRANSIT OBSERVATIONS 

5.1. Inner Planet: Transit Timing Scheduling 

As noted in section |4] from Fig |4] and Fig \5\ we have a potential means of distinguishing between different relative 
inclination and/or (outer) eccentricity states of the system. To do so will require a detailed investigation of the 
inner transit times to an accuracy greater than ~ 5 seconds. In this section we focus on the time-frame close to the 
pericenter- approach of the outer planet. It is at this point that the TTVs induced on the inner planet are predicted 
to go through a sharp change in values in only a small number of orbits. In Fig|5]we display the TTVs expected for 
three different system configurations: (i) irei = 0° and Couter = 0.69; (ii) irei = 45° and Couter = 0.69; (i) irei = 0° 
and Couter = 0.73, in which we have fixed the longitude of ascending node to be = 0, and chosen the other system 
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Fig. 6. — Contour plots in the {aouter /ciiriner , Pouter) plane, illustrating the contours of the predicted RMS TTV amplitude at different 
relative inclinations. The horizontal and vertical lines give the approximate observational constraints from Fig |3] (0.5 < Couter < 0.85, 
420 < Pouter < 440 days). The plots make clear that the TTVs can always be very large towards the top of the allowed eccentricity 
range, but that there is little variation with semi-major axis. Very high inclinations result in TTVs that are significantly greater than the 
observational constraints from B09 across the entirety of the parameter space. 



Jitter and TTVs in the HAT-P-13 system 15 




Relative Inclination (i^ei) 

Fig. 7. — Contour plots in the (irei, Pouter) plane, illustrating the contours of the predicted RMS TTV amplitude (The contours used 
are the same as those in Fig|6)l. The horizontal lines indicate the same upper and lower bounds as given in Fig. 6, as well as the best fit 
value for Couter from the right-hand column of Table 2, WIO. We can clearly see that, irrespective of the eccentricity of the outer planet, 
inclinations 88° < ij-ei < 92° will produce TTVs larger than the RMS constraints from Bakos ct al. (2009). Moreover, it is clear that if 
these observational constraints can be tightcncd-up, the allowed range of {ireh ^outer) will be significantly reduced, allowing them to be of 
significant more practical use. 
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Fig. 8. — Using transits of the inner planet to determine system inclination. We show examples of TTVs for the cases of ire.i = 0° and 
Vel = 45° favored by the Mardling (2010^) analysis. In the inset we provide TTV profiles over approximately 3 years, but in the main plot we 
focus on the region close to pericenter passage, illustrating that over the course of ~ 50 days (~ 17 transits), the variation chages from ~ —30 
to ~ 15. We suggest that sufficiently precise observations over this period should be able to distinguish the different inclination states. We 
also plot an example of an Vei = 0° with a different eccentricity (0.73), illustrating that sufficiently precise observations can also distinguish 
eccentricity states. We note that these illustrative simulations were chosen to have their time of pericenter passage at BJD = 2455302.2, 
the best fit date from Table [T] for 2.0 < o-j < 4.5 ms~^ and TTV amplitude less than 100 seconds: given the uncertainty in our fits for 
Tp^outer (see section [521l , the absolute dates may shift, but the general message regarding the relative timing of the large changes in signal 

amplitude will not. All plots have dinner — 0.043 AC/, Cinner — 0.02, Tninner — 0.85A^j, Ciouter — 1.188, &l TTlouter ^^^iouter = 15. 2M,/.) 



parameters to be those of the best fit value from the far right-hand column of Table 1. 

We see that in all three cases, the TTV signal changes sharply during the time between the predicted transit of the 
outer body (see section 15.21 below) and the predicted date of pericenter passage for the outer body. From Julian date 
^ 2455300, over the next ~ 50 days (i.e. over ^ 17 transits of the inner planet), the signal would be expected to switch 
from being ~ —30s to being ~ -|-15s. A similar pattern would be expected to occur at potential transit dates in 2011 
& 2012 at or around the dates given in Table 2. If enough observations can be made of high enough precision over this 
period, it should be possible to (i) At the very least confirm that a TTV variation is seen, (ii) give some significantly 
improved constraints on the magnitude of the TTVs in the system and hence point towards a better understanding of 
the system parameters, as well as (iii) pointing the way towards a subsequent, longer term observational campaign to 
further track the TTV profile over the course of an entire orbit of the outer planet. 

We would thus suggest that an observational campaign targeting the transits of the inner planet could give informa- 
tive results about the relative inclination of the planets in the system. We suggest that observations extending over 
the period of pericenter passage (for the outer planet) will initially be the most useful in determining the nature and 
amplitude of the TTVs, but that longer term observations can also give important information that allow different 
inclination states to be distinguished. In addition, it should be emphasized that such obse rvations would need to be 
of a precision s imilar to (or better than) that demonstrated in the study of WASP-10 bv iJohnson et all (|2009D and 
iJohnson et alJ ()2010( ). where the mid-transit times were determined to a precision of '-^ 7 seconds. If such a level of 
precision could be achieved over the course of 10 - 15 transits of the inner planet, then it is plausible that significant 
inclination information could be extracted. 

We caution that the precise date(s) of the pericenter passage for the outer planet are, of course, highly dependent on 
the values of the fitted orbital parameters aouter, Pouter, etc. Given that we have demonstrated in §3 that these fitted 
parameters are rather uncertain, it should be clear that the precise time of pericenter passage is similarly uncertain, 
and so any scheduled observations should take these uncertainties into account. Moreover, if/when further additional 
RV observations become available, it would be prudent to repeat the analysis of §3 to allow the uncertainties in the 
planetary parameters to be improved, and hence allow any future observational strategy to be refined. We also caution 
that in the examples shown in Fig [8l we have made specific assumptions regarding the longitude of ascending node of 
the planets (fixed = 0). If such assumptions are relaxed, then TTV profiles can be generated that differ in the degree 
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and extent to which their amphtude changes over the course of the pericenter passage. We emphasize that in fitting 
/ modeling any putative future observations, such degeneracies would need to be accounted for. 

5.2. Outer Planet: Transit Timing Projections 

It is clear from sections H] and 15.11 that further observations of the precise timings of the transit of the inner planet 
will be the key to understanding the TTVs, and subsequently the relative inclinations, of the planets around HAT- 
P-13. However, we feel that is is worth emphasizing that a (fortuitous) transit observation of the outer planet would 
help to pin down the inclination, period and most other parameters of the planetary system much more precisely (to 
say nothing of giving the first ground-based observation of a system with two transiting planets!). While the a-priori 
probability of HAT-P- 1 3c transiting is relativ e ly low if all inclination angles are considered equally likely (~ 1%, see 
iSeagroves etldl pOOl : E ane and von BraunI (|2009[ )) , this probability can b e significantly increased (to ^ 8%) for 
certain particular inclination and observation angles ()Beattv and Seageill201Cl[ ). 

If the outer planet does indeed transit, then there remains the issue of knowing when to observe the star to detect 
the transit. If we use our MCMC data to predict the times of the next few transits, assuming that the o uter p lanet is 
inclined at 90° to the plane of the sky (as was also assumed in B09), then using the method outlined in iKand ((2007); 
iKane and von BraunI ()2009[ ). we find that there is a significant uncertainty in the predicted timing of any potential 
transit. We note at this point that our potential transit date for the year 2009 {Tx^outer.Bakos, Table 1 & 2) is typically 
a few days later than that derived by B09 (but with a large range of uncertainties) and that our fitted outer period is 
slightly longer (see Table [Ij. 

Specific predictions require us to make some assumptions about the stellar jitter, so as an example we consider a low 
stellar jitter case and select the results in which the overall jitter lies in the range 2.0 < aj < 4.5 ms~^ (and take the 
cases with TTV amplitude less than 100 seconds). In this case, we find that the predicted transit dates for 2010 - 2012 
are TT,outer,20iQ = 2455314.3+?:^, TT,outer,20ii = 2455764.3^1^:^ and TT,outer,20i2 = 2456213.7tl^ j respectively, which 
translates into April 27th 2010 t^ j days, July 21st 2011 tlH days and Oct 13th 2012 tlfl days. We display in Fig[9] 
the scatter in the predicted transit times for 2011 as a function of the assumed jitter and an accompanying histogram 
to illustrate the spread in predicted transit times, as well as the difference that results from taking a restricted range 
of jitter values. We do this for both the B09 SI and WIO S2 data sets. It is clear from the scatter plot that within 
a given data set (e.g. B09) the distribution is essentially the same as that of the outer period plotted in Fig [T] The 
histograms suggest that the systems with jitter in the range 2.0 < aj ^ 4.5 m s~^ tend to have predicted transits which 
are more sharply peaked towards slightly lower values than the distribution which includes all jitter values. 

However, when one compares the SI results to those of S2, we see that there is an offset in the predictions, due to 
the preference for longer periods for the outer planet in S2, as well as the asymmetry in the parameter space. 

We again emphasize that different assumptions regarding the level of jitter in the system act to shift the predicted 
date (as well as the error bars) of any potential transit. We also stress that the uncertainties on our figures cover 
only two-thirds of the data: more extreme excursions exist. It would thus seem prudent for any future campaign of 
transit observations to be undertaken significantly before and after the mid-point of any predicted transit date, to 
allow for the uncertainties: e.g. for the 2011 potential transit, don't just observe on July 21st, 2011, but rather take 
observations over a large range of nights (perhaps starting around July 9th and continuing for ~ 3 weeks) to give the 
highest probability of observing any transit that might occur. 

On a practical point, we note from the object visibility calculator at ' http:/ /catserver. ing.iac.es/staralt/ ' that HAT- 
P-13 will be visible for many hours per night in the northern hemisphere during October 2012, but that observations 
around July 2011 would be rather more challenging, as the star is likely to only be visible for ^ 1 hour per night. 

5.3. Additional Thoughts on Variations in Cinner o-nd rj — vjinner — Pouter 

We note from the study by iMardlingI (|2010( ) of the likely dynamics of the HAT-P-13 system, that the secular 
timescale for Cmner and rj to vary around the limit cycle is of the order a few thousand years (see her Fig. 3b). We 
can also see that the amplitudes of the oscillations are ^ 50° for rj and ^ 5 x lO"'^ for einner- This would suggest that 
over the course of a 10 year observational campaign, the change in these quantities would be Ar/(10j/rs) ^ 0.2° and 
Aeiri„er(10yrs) ~ 10~^, both very small quantities indeed. 

6. CONCLUSION 

We have performed an investigation of the HAT-P-13 system, concentrating our efforts on re-evaluating the RV 
analysis, including jitter as an intrinsic part of the MCMC analysis, and then combining this evaluation with an 
analysis of the transit timing variations (TTVs) in the system. 

We find that: 

• If the methodology for including jitter within an RV analysis is changed and jitter is included as a model 
parameter within the MCMC analysis, a significantly larger range of parameter space opens up, i.e. the masses 
and orbital parameters for the system are significantly less well-defined. As an example, based on the extended 
RV data set analyzed in Winn et al. 2010, if the overall system jitter is in the range 2.0 < aj < 4.5 m s~^ rather 
than aj ~ 3.4 ms""'^, the eccentricity of the inner planet has a one-sigma best fit value of emner = 0.038lQQ^g, 
the eccentricity of the outer planet has best fit value of 0.65to;j!|g, while the period of the outer planet would have 
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Fig. 9. — Predicted transit dates for HAT-P-13c using the MCMC sample fits from B09 (top) and WIO (middle). We assume an inclination 
for the orbit of 90° to the plane of the sky and restrict our attention here to only those systems which have TTVs less than 100 seconds. 
The dates are Julian dates -2450000. Top & Middle It is clear that the transit time distribution is of the same form as that of the period 
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different jitter levels are similar within a data set, but that the different data sets can have the center of the distributions shifted relative 
to one another due to the preference for a larger Pouter in the S2 analysis 
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best fit values of Pouter = 449.4^42 and the relative pericenter alignment of the two planets becomes essentially 
unconstrained, with a best fit value of 77 = 12.4^g3g. With the current data set, even higher jitter values are 
plausible, and this would act to further increase the uncertainty in the observations of system parameters (see 
Section EH). 

• If we include the current weak TTV constraints (< 100 seconds) in addition to the RV analysis, then we find that 
we can exclude eccentricities for the outer planet larger than ^ 0.85. If future observations can pin-down the 
RMS TTV amplitude to a much narrower range (< 10 seconds for example), then this would strongly constrain 
the eccentricity of the outer planet ( < 0.6 for the case of TTVs less than 10 seconds.). See Section l3^ for a 
discussion of the implications of the TTV amplitude. 

• The current TTV constraints already suggest that the planets in the system do not have relative inclinations in 
the range 88° < irei < 92°, but any future tightening of the TTV constraint would act to exclude a much wider 
range of inclinations centered on i^ei = 90° (70° < Vef < 110° would be excluded if TTVs amplitudes were less 
than 10 seconds - See Section |4?2|) 

• The TTV profile can in many circumstances act as an efficient diagnostic tool in determining the relative 
inclination between the two planets. In particular, it is easy to distinguish be tween s ystems with i ^ 0° and 
i 45° (the approximate angles found to be the most likely in the analysis of Mardlingl () 20101) ) through suitably 
timed transit observations (see Section [4.1^ . We thus suggest that an observational campaign targeting the 
transits of the inner planet could give informative results about the relative inclination of the planets in the 
system. We suggest that observations extending over the period of pericenter passage (for the outer planet) 
will initially be the most useful in determining the nature and amplitude of the TTVs, but that longer term 
observations can also give important information that allow different inclination states to be distinguished. 

• We note that a transit by the outer planet is rendered more likely following the exclusion of regions of inclination 
space by the analysis of iMardling (|2010l ) , supporting the idea of an observational campaign, but the large 
spread in allowed orbital periods resulting from our jitter analysis would suggest that any observations should 
be conducted over a ~ 3 week period centered approximately on either July 21st 2011, or October 13th 2012. 
(Observations in 2011 will be significantly more difficult than in 2012 - See Section [5]) 

• Given the value of transit timing observations near the time of pericenter of the outer planet, the long orbital 
period of the outer planet, and the poor observability during pericenter passages in 2011, we suggest that 
observers should consider observing whenever it would be possible to observe full ingress or full egress, rather 
than only observing full transits. Measuring a large fraction of transits/ingress/egress times during this critical 
time will require observations from a network of telescopes at multiple longitudes. 
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Fig. 10. — TTV Amplitude vs. System Parameters for parameters in which there is no obvious correlation between amplitude and 
parameter value. In the various plots above we show the RMS TTV amplitude as a function of various system parameters (planetary 
eccentricity, etc). As is Fig. 1, systems with Cj < 2.0ms~^ are plotted using a red circle (gray in the print version), those with 
2.0 < (Tj < 4.5ms~^ are plotted using a green cross (gray in the print version), and those with aj > 4.5ms~^ are plotted using black 
addition symbols. As in Fig. 1, in the left-hand column we present results using only the subset of data known at the time of publication 
by B09, assuming a 2-planet fit. In the center and right-hand columns we present results obtained using the full data set of WIO, with the 
central column assuming 2-planets + a linear trend, while the right-hand column assumes a 3-planet fit. We plot results for the period of 
the outer planet. Unlike the results plotted in Fig[2l we see that there is no obvious correlation between the RMS TTV amplitude and the 
system parameter value. 
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APPENDIX 

For the sake of completeness, we return to consider Fig [5] in which we demonstrated that the fitted values for the 
overall RV offset, as well as the RV Amplitude, K outer , longitude of pericenter, w outer, and the eccentricity, eouter, of 
the outer planet all showed a correlation with the RMS TTV amplitude. In Fig. 10 we plot Pouter for which there is 
no obvious correlation between the RMS TTV amplitude and the fitted system parameter. 



